Inferring coupling strength from event-related dynamics 
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We propose an approach for inferring strength of coupling between two systems from their tran- 
sient dynamics. This is of vital importance in cases where most information is carried by the 
transients, for instance in evoked potentials measured commonly in electrophysiology. We show 
viability of our approach using nonlinear and linear measures of synchronization on a population 
model of thalamocortical loop and on a system of two coupled Rossler-type oscillators in non-chaotic 
regime. 
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I. INTRODUCTION 

Coherent actions of apparently distinct physical sys- 
tems often provoke questions of their possible interac- 
tions. Such coherence in interacting systems is often a 
result of their synchronization [ll]. It became a popu- 
lar topic with the discovery of synchronization of non- 
identical chaotic oscillators Over the years dif- 
ferent types of synchrony were studied, notably phase 
synchronization fsl]. There were also numerous at- 
tempts to study more compHcated interactions under the 
names of generalized synchronization or interdependence 
[1, H, @, 0, Hi- In biological context synchronization is 
expected to play a major role in cognitive processes in 
the brain (9l.fl0l.Tll| such as visual binding and large- 
scale integration [ll|. Various synchronization measures 
were successfully applied to electrophysiological signals 
[ul, 113, H El, m, lli, [13, El- In this work we concen- 
trate on nonlinear interdependence [l3. fl3|. 

For an experimentalist it is often interesting to know 
how two systems synchronize during short periods of 
evoked activity [H, [l^l ■ Such questions arise naturally in 
analysing data from animal experiments (ill, lil. [2^. E^l . 
One measures there electrical activity on different levels 
of sensory information processing and aims at relating 
changes in synchrony to the behavioral contex, such as at- 
tention or arousal. It may be the case that the stationary 
dynamics (with no sensory stimulation) corresponds to a 
fixed point. For instance, when one measures the activity 
in the barrel cortex of a restrained and habituated rat, 
the recorded signals seem to be noise (2ll. |2^. |23| . On the 
other hand transient activity evoked by specific stimuli 
seems to provide useful information. For example, bend- 
ing a bunch of whiskers triggers non-trivial patterns of 
activity (evoked potentials, EPs) in both the somatosen- 
sory thalamic nuclei and the barrel cortex (23l . [25| . 

Explorations described in this paper aim at solving the 
following problem. Suppose we have two pairs of tran- 



sient signals, for example recordings of evoked potentials 
from thalamus and cerebral cortex in two behavioral situ- 
ations [2ll.l23|. Can we tell in which of the two situations 
the strength of coupHng between the structures is higher? 
Thus we investigate if one can measure differences in the 
strength of coupling between two structures using non- 
linear interdependence measures on an ensemble of EPs. 
Since EPs are short, transient signals, straightforward 
application of the measures motivated by studies of sys- 
tems moving on the attractors (stationary dynamics) is 
rather doubtful and a more sophisticated treatment is 
needed [13, HH]- Our approach is similar in spirit to that 
advocated by Janosi and Tel for the reconstruction of 
chaotic saddles from transient time series [23| • (Note that 
the transients we study should not be confused with the 
transient chaos studied by Janosi and Tel.) Thus we cut 
pieces of the recordings corresponding to well-localized 
EPs and paste them together one after another. Since 
we are interested in the coupled systems, unlike Janosi 
and Tel, we obtain two artificial time-series to which we 
then apply nonlinear interdependence measures and lin- 
ear correlations. It turns out that this approach allows 
to extract the information about the strength of the cou- 
pling between the two systems. 

We test our method on a population model of informa- 
tion processing in thalamocortical loop (Figure [T|) con- 
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Figure 1: Structure of the model of the thalamocortical loop 
used in the simulations. 
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sisting of two coupled Wilson-Cowan structures |28l. l29l|. 
Sensory information is relayed through thalamic nuclei 
to cortical fields, which in return send feedback connec- 
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tions to the thalamus. This basic framework of the early 
stages of sensory systems is to a lar ge e xtent universal 
across different species and modalities [30|. To check that 
the results are not specific to this particular system we 
also study evoked dynamics of two coupled Rossler-type 
oscillators in non-chaotic regime. 

The paper is organized as follows. In Sec. Ull we de- 
fine the measures to be used. In Sec. [Till we describe the 
models used to test our method. Our model of thalamo- 
cortical loop is discussed in Sec. lIII Al and a system of two 
coupled Rossler-type oscillators is described in Sec. lIIIBl 
In Sec. HVlwe present the results. In Sec. IIV Al we show 
how various interdependence measures calculated on the 
transients are related to the coupling between the sys- 
tems, while in Sec. IIVBI we study how the resolution of 
our methods degrades with noise. Finally, in Sec. IIV C[ 
we apply time-resolved interdependence measure Hi [2a| 
and compare its utility with our approach. We summa- 
rize our observations in Sec. |Vl 



temporal correlations [3J|, points closer in time to the 
current point x„ than a certain threshold are typically 
excluded from the nearest-neighbor search (Theiler cor- 
rection). Then we define the y-conditioned mean 

i?«(X|Y) = i^(x„-x,„j2^ 



where the indices r„ j of the nearest neighbors of x„ are 
replaced with the indices s„ j of the nearest neighbors of 

y„. The definitions of Rn\Y) and i?i'^''(Y|X) are anal- 
ogous. The measures H and N use the mean squared 
distance to random points: 



i?„(X) = 



and are defined as 



N 



II. SYNCHRONIZATION MEASURES 

In the present paper we mainly study the applicability 
of nonlinear interdependence measures on the transients. 
These measures, proposed in are non-symmetric and 
therefore can provide information about the direction of 
driving, even if the interpretation in terms of causal re- 
lations is not straightforward [sij. 

These measures are constructed as follows. We start 
with two time series Xn and ?/„, n = 1, . . . , iV, measured 
in systems X and Y. We then construct m-dimensional 
delay-vector embeddings [s^l x„ = {xn, ■ ■ ■ , a;„_(m_i)T-), 
similarly for y„, where r is the time lag. The informa- 
tion about the synchrony is inferred from comparing the 
size of a neighborhood of a point in m-dimensional space 
in one subsystem to the spread of its equal-time coun- 
terpart in the other subsystem. The idea behind it is 
that if the systems are highly interdependent then the 
partners of close neighbors in one system should be close 
in the other system. Several different measures explor- 
ing this idea can be considered depending on how one 
measures the size of the neighborhood. These variants 
include measures denoted hy S, H N M 
We have studied the properties of most of these measures 
but for the sake of clarity here we report only the results 
for the "robust" variant H and a normalized measure N, 
as they proved most useful for our purposes. 

Let us, following [Tl|, for each x„ define a measure of 
the spread of its neighborhood equal to the mean squared 
Euclidean distance: 



where r„ are the time indices of the k nearest neighbors 
of x„, analogously, Snj denotes the time indices of the 
k nearest neighbors of y„ . To avoid problems related to 
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The interdependencies in the other direction i/'^'^^ (Y|X), 
7V('=)(Y|X) are defined analogously and need not be equal 
i/W(X|Y), iVW(X|Y). 

Such measures base on repetitiveness of the dynam- 
ics: one expects that if the system moves on the attrac- 
tor the observed trajectory visits neigborhoods of every 
point many times given sufficiently long recording. The 
same holds for the reconstructed dynamics. However, if 
the stationary part of the signal is short or missing, espe- 
cially if we observe a transient such as evoked potential, 
this is not the case. Still, if we have noisy dynamics, every 
repetition of the experiment leads to a sHghtly different 
probing of the neighborhood of the noise-free trajectory. 
This observation led us to an idea of gluing a number 
of repetitions of the same evoked activity (with differ- 
ent noise reaHzations) together and using such pseudo- 
periodic signals as we would use trajectories on a chaotic 
attractor. A similar idea was used by JanjDsi and Tel in 
a different context for a different purpose ||27]. An exam- 
ple of a delay embedding of a signal obtained this way 
is presented in Fig. [2l Note that artifacts may emerge 
at the gluing points. This is discussed in [131 , and some 
countermeasures are proposed. For simplicity we proceed 
with just gluing as we expect that the artifacts only in- 
crease the effective noise level. The infiuence of noise is 
studied in Sec. IIVBI 

Recently, time-resolved variants of the methods de- 
scribed above were studied (20l . |2^ . They are applied 
to ensembles of simultaneous recordings, each consisting 
of many different realizations of the same (presumably 
short) process. Let us denote the n-th state vector in j- 
th realization of the time-series by X:^ (y^ , respectively) , 
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material to [26|. To compare the linear and nonlinear 
analysis methods we calculated the cross-correlation 
coefficients using Matlab. 

While in numerical studies the correctness of recon- 
struction can often be easily checked by comparison with 
original dynamics, in analysis of experimental data it can 
be a complex issue. Correct reconstruction is a prerequi- 
site for application of our technique. For technical details 
on best practices of delay embedding reconstructions, pit- 
falls and caveats, see |35l |. 



Figure 2: Delay-vector embeddings (shown in planes defined 
by the first two principal components) of pseudo-periodic sig- 
nals obtained by gluing 50 evoked potentials generated in a 
model of thalamocortical loop. On the left (signal from "tha- 
lamus") a point is chosen (black square) and its 15 nearest 
neighbors are marked with red (gray) diamonds. On the right 
("cortex") the equal-time partners of the marked points from 
the left picture are shown. 



j = 1, . . . , J. The idea in ^] is, for given to find 
one neighbor in each of the ensembles. Then a mea- 
sure (denoted T) based on distances to these neighbors 
is constructed. The proposition of [2^1 is to look not at 
the nearest neighbors of a given x„ no matter what time 
they occur at, but rather at the spread of state- vectors 
at the same latency across the ensemble. In Sec. IIV CI 
we study the measure Hi as defined in jH]. Let r-'' de- 
note the ensemble index of the ^-th nearest neighboor of 
y:^ among the whole ensemble {y^i}-'^^' ' '^. Define the 
quantities 

i?fW(X|Y) = i^(xi-x:'-v, 



R: 



III. MODEL DATA 

A. Connected Wilson-Cowan aggregates 

Our model of the thalamocortical loop was based on 
the Wilson and Cowan mean-field description of inter- 
actin g p opulations of excitatory and inhibitory neural 
cells [23, l29|. In the simplest version, which we used, 
each population is described by a single variable stand- 
ing for its mean level of activity 



dE 

'^E-^ = -E+ (kE - rEE)SE{cEEE - cjeI + P), 

Ti^ = -I+ {ki - riI)Si{cEiE - ciil + Q). 
at 



(1) 



The time-resolved interdependence measure is further de- 
fined as 

gW(X|Y) = lVlog . 

Analogously one can define i/|'^-'(Y|X) and also time- 
resolved variants of other interdependence measures. 

In the numerical experiments described in this paper 
we use the following parameters of the nonlinear interde- 
pendence measures: time lag for construction of delay- 
vectors: T = 1, embedding dimension m = 10, number 
of nearest neighbors A: = 15, Theiler correction T = 5. 
To calculate the interdependencies we used the code by 
Rodrigo Quian Quiroga and Chee Seng Koh available at 
|http : / / www . vis . caltech . edu/~rodri/Synchro/Synchro, 
In case of the measure Hi we use the same embedding 
dimension and time lag; here k = \. To calculate this 
measure we used the code provided in supplementary 



The variables E and / are the mean activities of excita- 
tory and inhibitory populations, respectively, and form 
the phase space of a localized neuronal aggregate. The 
symbols r, fc, r, c denote parameters of the model, S 
are sigmoidal functions, P and Q are input signals to ex- 
citatory and inhibitory populations, respectively. These 
equations take into account the absolute refractory pe- 
riod of neurons which is a short period after activation 
in which a cell cannot be activated again. Such models 
exhibit a number of different behaviors (stable points, 
hysteresis, limit cycles) depending on the exact choice 
To relate the simulation results 
23l we considered the observable 



of parameters 2^, 
to the experiment 
V — E — I , since the electric potential measured in ex- 
periments is related to the difference between excitatory 
and inhibitory postsynaptic potentials (see the discussion 
in (H). 

We studied a model composed of two such mutually 
connected aggregates, which we call "thalamus" and "cor- 
tex" (Figure [Ij. Note that the parameters characterizing 
the two parts are different (see the Appendix □ for a com- 
plete specification of the model). Specifically, there are 
no excitatory-excitatory nor inhibitory-inhibitory con- 
nections in the thalamus. Only the thalamus receives 
sensory input, and we assume that Q is always a constant 
fraction of P. The connections between two subsystems 
are excitatory only. 

.hoflfe. Moiiel the stimulus we assumed that the input 
(P, Q) switches at some point from to a constant 
value {PcQc)^ and after a short time (on the time- 
scale of relaxation to the fixed point) switches back to 



4 



zero. This is clearly another simplification, as the real 
input, which could be induced by bending a bunch of 
whiskers [U, [H, [IS] , would be a more complex function 
of time. However, the transient nature of the stimulus 
is preserved. In this simple setting we can understand 
that the "evoked potential" corresponds to a trajectory 
approaching the asymptotic solution of the "excited" sys- 
tem (with the non-zero input Pc, Qc), followed by a re- 
laxation to the "spontaneous activity" in the system with 
null input. 

The model parameters were chosen so that its response 
to brief stimulation were damped oscillations of V both 
in the thalamus and the cortex similar to those observed 
in the experiments, both in terms of shape and time du- 
ration [21I, m, [131 (Figure [3]). However, apart from that, 
we exercised little effort to match the response of the 
model to the actual activity of somatosensory tract in 
the rat brain. Our main goal in the present work was es- 
tablishing a method of inferring coupling strength from 
transients and not a study of the rat somatosensory sys- 
tem. For this reason it was convenient to use a very sim- 
plified, qualitative model. Interestingly, the response of 
the model, measured for example as the activity of exci- 
tatory cells in the thalamus, extends in time well beyond 
the end of the stimulation (Figure [3]) . Such behavior is 
not observed in a single aggregate and requires at least 
two interconnected structures [29| . 



(a) (b) 
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lO^V^^(PC^) 10'V^^(PC,) 

Figure 3: "Evoked potentials" (V = E - I), (a), (b) and 
their delay-vector embeddings shown in a plane defined by 
the first two principal components (c), (d). Plots (a) and (c): 
thalamus, (b) and (d): cortex. The intervals above the EP 
indicate the duration of the non-zero stimulus. Black (thick) 
lines are solutions for the system without noise, blue (thin) 
curves are five different realisations of noisy dynamics. 

We performed numerical simulations in three modes: 
either stationary (null or constant input), or not (tran- 
sient input). The dynamics of the model is presented in 
Figure [4l In case of transient input the simulation was 



done for -1000 < t < 1000ms. We used the stimulus P 
and Q which was except for the time 200 < t < 220 
when it was Pc = 3.5 and Qc — 0.3. The system set- 
tled in the stationary state during the initial segment 
{t < 195) which was discarded from the analysis. The 
noise was simulated as additional input to each of the 
four populations, see the Appendix □ for the equations. 
For each population we used different Gaussian (mean 
fi — 0, standard deviation a = 0.025) white noise, sam- 
pled at IkHz and interpolated Hnearly to obtain values for 
intermediate time points. In case of stationary dynamics 
we simulated longer periods, —1000 <t< 20000ms. The 
signals were sampled at lOOHz before the synchronization 
measures were applied. 

In case of constant or null stimulation the system ap- 
proaches one of the two fixed-point solutions which are 
marked by large dots in Figure |4l For the amount of 
noise used here the dynamics of the system changes as 
expected: the fixed points become diffused clouds (Fig- 
ure |4]) . During the transient — "evoked potential" — 
the switching input forces the system to leave the null- 
input fixed point, approach the constant-input attractor, 
and then relax back to its original state (Figure [4]) . Of 
course, in the presence of noise the shape of the transient 
is affected (Figured]). Observe the similarity between the 
embedding reconstructions of the evoked potentials (Fig- 
ure [3l bottom row) and the actual behavior in Vrh-^cx 
coordinates (Figure [4l bottom row). 



B. Coupled Rossler-type oscillators 

While we are specifically interested in the dynamics 
of thalamocortical loop which dictated our choice of the 
studied system, we checked if our approach is not specific 
to this model. Our second model of choice consisted of 
two coupled Rossler-type oscillators 0, Is^ 

—— = -(1 + Aa;)i/i - zi + aC(a;2 - +<^i, 
at 

% = (1 + Acj)xi -0.15yi + P + 6, 
at 

^ = 0.2 + zi(xi -10) + C3, 
at 

—— = -{1 - Auj)y2 ~ Z2 + aC{xi ~ X2) + 
at 

% = il-Auj)x2-0.15y2 + ^5, 
at 

= 0.2 + Z2(x2-10)+C6. 

We used the frequency detuning parameter Aoj = 0.05 
and the maximum coupHng constant C = 0.06. The scal- 
ing parameter a took values from to 1. The stimula- 
tion parameter P was except for 200 < i < 250 where 
it was set to 0.8; the noise inputs £,i, i = 1...6 were 
Gaussian white noise with parameters as for the Wilson- 
Cowan model. The simulation was done for t g [0, 300] 
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Figure 4: Dynamics of the model. The green (lower left) and red (upper right) dots are fixed points in case of null or constant 
stimulation respectively, the black (thick) line is the noise-free transient dynamics. Blue (thin) lines are example trajectories 
of the model in the presence of noise. The plots show projections of the same dynamics to different planes. 



and segments from t = 195 to t ~ 300, sampled every 
At = 0.125, were used for the analysis of the transients. 
The synchronization was measured between xi and X2. 
Parameters of the system were chosen so that asymptot- 
ically it moved into a stable fixed point (note the signs in 
the equations for yi and ^2) for both P = and P = 0.8. 
Therefore the transient dynamics (Fig. [5]) is of the same 



(a) (b) 




X, (PC^) XjlPC,) 



Figure 5: (a), (b): signals (a; coordinate) from coupled 
Rossler-type oscillators; (c), (d): their delay- vector embed- 
dings, shown in a plane defined by the first two principal 
components. The intervals in (a) and (b) indicate the dura- 
tion of non-zero input P. Black (thick) lines are solutions for 
the system without noise, blue (thin) curves are five different 
realisations of noisy dynamics. 

type as in the model of thalamocortical loop: the sys- 
tem switches briefly to the second stable point and then 
returns. Note that the level of noise in the second subsys- 
tem is quite high and the evoked activity is barely visible 
at the single trial level (Fig.[5l right column). 



IV. RESULTS 

A. Inferring connection strength 

We aim at solving the following problem: suppose we 
have two pairs of signals, for example recordings from 
thalamus and cerebral cortex in two behavioral situa- 
tions [m, m, [13, HI- Can we tell in which of the two 
situations the strength of connections between the struc- 
tures is higher? Thus we need to flnd a measure being 
a monotonic function of the coupling strength. We have 
studied this problem in our model of thalamocortical loop 
(Section IlII Ap . We scaled the strength of connections 
from thalamus to cortex by changing a between and 1, 
and calculated synchrony measures on signals from these 
structures. The strength of connections from cortex to 
thalamus was constant {(3 = 1); see the Appendix □ for 
the details. 

Consider first stationary signals with P = or P = 
const. Without noise the system is in a fixed point 
and obviously it is impossible to obtain the connection 
strength. However, given the noise, in principle the dy- 
namics in the neighborhood of the fixed point is also 
probed. Thus there is a possibility that the interdepen- 
dence and the strength of the coupling could be estab- 
lished during stationary parts of the dynamics. It turns 
out that for null stimulation neither the interdependence 
measures nor the linear correlations detect any changes in 
the coupling strength (Figure [6l left column) . For con- 
stant non-zero input there is a connection between the 
coupHng strength and the values of the measure but they 
are anti-correlated and the dependence is not very pro- 
nounced (Figure [6l right column). One must also bear in 
mind that while it is possible to have no stimulation, in 
brain studies prolonged and constant stimulation in the 
present sense cannot be experimentally realized (at least 
for most sensory systems) because of the adaptation of 
receptors. The natural stimuli are necessarily transient. 

To use the synchrony measures on the transient we 
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Figure 6: Nonlinear interdependence measures (a) H, (b) A'^ and maximum of absolute value of cross-correlation coefficients 
for varying coupling constants. Each measure is calculated for 10 independent realizations with different seeds. The parameter 
a scales the strength of connections from thalamus to cortex. 



cut out pieces of signal corresponding to the evoked po- 
tential, and pasted them one after another. Thus ob- 
tained pseudo-periodic signal contained the same under- 
lying dynamics with each piece differing due to the noise. 
We then appHed the same measures as we did for the 
stationary signals. In the simulations we calculated 50 
"evoked potentials" (Figure [3]) for each value of a. Plots 
in the middle column of Figure [6] show the values of the 
synchronization measures evaluated for different coupling 
strengths. It can be seen that they are increasing func- 
tions of the coupling strength between the subsystems. 
Therefore, our approach is indeed a viable solution to 
the problem of data-based quantification of the coupHng 
strength. 

It is interesting to study the values of these inter- 
dependence measures in different cases. Observe that 
H{VTh\Vco,) > i?(Vcx|VTh) for P = 0. The opposite 
is true for transients (for small a). This is even more 
clearly visible for TV. In all the cases linear correla- 
tions showed similar trends to the nonlinear measures 
N{VTi,\Vc^),N{Vc.\VTh)- 

The asymmetry in the interdependence measures was 
originally intended to be used for inferring the direction 
of the coupling or driving. However, the inference of spe- 



cific driving structure in every case must follow a careful 
analysis of underlying dynamics (see, for example, dis- 
cussions in [sJl and Let us consider the plots in 
the middle column of Figure [6l For small a the domi- 
nant connections are from the cortex to the thalamus so 
one might expect that the state of the thalamus might 
be easier predictable from the states of the cortex than 
the other way round. Thus one would intuitively expect 
i/(VTh|Vcx) > H{Vcx\VTh)- However, we observe the 
opposite. The reason is that the measures used are re- 
lated to the relative number of degrees of freedom [l^ . 
Loosely speaking, as discussed fsij, the effective dimen- 
sion of the driven system (thalamus for small a) is usu- 
ally higher than the dimension of the driver (which means 
that the response — the dynamics of the thalamus — is 
"more complex"). This effect is further enhanced by the 
fact that we stimulate the thalamus in moments unpre- 
dictable from the point of view of the cortex. Summa- 
rizing, the result is compatible with the analysis in [3l] |. 
What happens for higher a when the two measures be- 
come equal is probably the coupHng between the two sub- 
systems becoming so strong that the quality of prediction 
in any direction is comparable. 

In the stationary case the situation is different as we 
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observe the asymptotic behavior. It turns out that for 
P = for every a, and for P = const > for small 
a we have H{VYh\Vcx) > ^^(^Cx|VTh)- But it seems 
that another effect also plays a role here. The noise in 
the cortex has a higher amplitude than in the thalamus 
and as a consequence it is easier to predict the state of 
the thalamus from that of the cortex than in the other 
direction. The reason for this disparity in the amplitudes 
is the difference in the shape of the sigmoidal functions 
Sq. To summarize, here, the asymmetry of the measures 
reflects internal properties of the two subsystems and not 
the symmetry properties of the coupling between them. 

Figure [7| shows similar results obtained for two coupled 
Rossler-type systems. In stationary situation the inter- 
dependence measures are very noisy. Although a weak 
trend is visible, one would not be able to reliably dis- 
criminate between, say, a = 0.25 and a = 0.75. The 
equality of the measures in two directions is due to the 
fact that the systems are almost identical and symmetri- 
cally coupled. 

If the interdependence is quantified on transient parts 
of the dynamics, the situation improves considerably. 
H{X2\Xi) has a high slope and is a very good measure 
of the coupling strength between the systems. Although 
i/(Xi|X2) has a slope comparable to that in the station- 
ary case for P = 0, the variability of the results is much 
smaller, compared to the size of the fluctuation in the 
ensemble mean in the stationary case. The difference be- 
tween i/(X2|Xi) and i/(Xi|X2) reflects the asymmetry 
of the driving (which makes the dynamics of Xi "more 
complex" than the dynamics of X2), not of the coupling 
(which is symmetric). 



B. Influence of noise 

The performance of the procedure described above de- 
pends on the level of noise present in the system. To 
study this dependence we performed the simulations of 
the thalamocortical model (the case of transient dynam- 
ics) for 25%, 50%, 100% and 200% of the original noise 
level. We found that for increasing level of noise the dy- 
namics of the system may change qualitatively: if the 
noise level is large enough the system may be kicked out 
of the basin of attraction of the flxed point and would 
not return there after P is reset to 0. Instead it may 
fall into the basin of attraction of another stable orbit 
or switch between the basins repeatedly. We observed 
such behavior only once for 2500 simulations performed 
with 200% of the original noise and this trial was ex- 
cluded from the analysis. Such behavior becomes more 
frequent with increasing noise (e.g. 400%) and so we did 
not study this situation as it was very different from the 
original dynamics of the system. 

As one would expect, the higher the noise, the less sen- 
sitive the measures are (Fig. [HI). However, even for twice 
the original level of noise a weak trend in the interdepen- 
dence is clearly visible. 



C. Time-resolved measure Hi 

Since we are interested in the dynamics of non- 
autonomous systems one might wonder if time-resolved 
measures, such as Hi introduced in pB], would not 
perform better in the problem of inferring connection 
strength. We performed tests on cut-and-pasted tran- 
sient signals. This problem is different from the one stud- 
ied in There, two Lorenz systems were coupled for 
short periods of time and Hi was shown to identify these 
times of coupling well. In our problem the coupling is 
constant in time, it is only the input to the system that is 
varying. For the problem at hand the values of Hi do not 
seem to change with varying coupling constant a (Fig. [9l 
(a)) when (3 is constant, = 1. The reason for this 
may be that even for a = the subsystems are coupled 
through the connections from cortex to thalamus. This 
hypothesis can be tested in another experiment, where all 
the connections between the subsystems are scaled and 
a = p. Indeed, in this setup the measure Hi is sensitive 
to the coupling strength (Fig. [9l (b); Fig. [T0|) . 

One may also note that Hi(VTh|Vcx) is on average 
higher than Hi{Vc^\VY-h), exactly as for H in case of 
P = and contrary to what is observed using H on tran- 
sients (Fig. El (a)). Thus it seems that for the problem of 
inferring coupling strength between two systems the op- 
timal approach is to use H or N, or linear correlations, 
on the transients, as described in Section [IV Al 



V. CONCLUSIONS 

To summarize, we have proposed a general approach 
for inference of the coupling strength using transient 
parts of dynamics. We have shown that our approach 
gives more information about the coupling between sub- 
systems than the approach using the stationary part of 
dynamics in case when the asymptotic dynamics is on a 
flxed point. We have checked the validity of this approach 
on a model of a thalamocortical loop of sensory systems 
and on two coupled Rossler-type oscillators. We showed 
that our method is quite robust with respect to increas- 
ing level of noise as long as the dynamics does not change 
qualitatively. We have also shown that this method mea- 
sures different aspects of coupling than a time-resolved 
measure Hi and than linear correlations. We believe that 
our approach will be of use in many other physical sys- 
tems studied in the stimulus-response paradigm, espe- 
cially in the experimental context. 

The results of Section IIV Al are compatible with our 
preliminary studies of data from real neurophysiologi- 
cal experiments 'j^^. There one cannot discern coupling 
strength in two contextual situations basing on station- 
ary recordings, but the analysis of transients leads to 
clear differences between two variants of experiment. The 
results of this analysis will be published elsewhere. 



8 



(a) Stationary, P=0 (b) Stationary, P=0, means (c) Transient 




Figure 7: Nonlinear interdependence measures H{Xi\X2), -H'(X2|Xi) and maximum of absolute value of cross-correlation 
coefficients between signals from two symmetrically coupled Rossler-type systems with noise. Coupling strength is proportional 
to a. The panels (a) and (c) present results of 10 simulations with different seeds. Additionally in stationary situation the means 
across repetitions are plotted for clarity in (b), the top curve (A) is cross-correlation. In stationary situation both nonlinear 
measures (B), (C) take similar values. On transients (c) H{X2\Xi) (C) is higher than H{Xi\X2) (B). The intermediate curves 
(A) are cross-correlation. 




0.5 



Figure 8: Nonlinear interdependence H in thalamocortical 
loop model for varying level of noise. The curves represent 
the values of the measure on transients. Noise level are 25% 
(a), 50% (b), 100% (c) and 200% (d) of the original noise. 



(a) 



a= 1, p = 1 



a = 0.27, p = 1 



a= 0, p = 1 



'4^;^ ' t^^^ 



500 1000 500 1000 500 1000 
Time [ms] Time [ms] Time [ms] 



(b) 



a = p = 1 



a=p = 0.27 



a = p = 



500 1000 500 1000 500 1000 
Time [ms] Time [ms] Time [ms] 



Figure 9: Time-resolved nonlinear interdependence measure 
^^^i(VTh|Vcx) (dashed lines) and Hi{Vcx\VTh) (solid lines) for 
3 different values of the coupling a. The means are shown 
with horizontal lines, (a) P — 1, (h) (3 — a. 
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We use the following equations for the model of thala- 
mocortical loop: 

xSi^,^ {Q + c2Et:i, + (3e2Ec^ + ^2), 
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1.5, 
1 

0.5 



0.5 




where 



Figure 10: Mean value of the time-resolved interdependence 
measure -ffi(VTh|Vcx) (dashed lines) and -ffi(Vcx|VTh) (solid 
lines), (a) = 1, (b) /3 = a 



dt 



dJ, 



Cx 



dt 



-Ecx + {kEc^ - rE'cx) 

--fcx + (^7cx - ^-^Cx) 

x5/c,(c5i^Cx - ce/cx + aaETi, + ^4), 



Sq{x) 



q standing for , -^Th , £^Cx , -^Cx , and ^i, i = 1 . . . 4 are 
noise inputs. The normalizing constants kg are defined 
as kg ~ 1 ~ j^^ga,fl, • 

In the numerical experiments we used the following 
parameter values: 



Cl 


= 1.35 


C2 


= 5.35 


C3 


= 15 


C4 


= 15 


C5 


= 15 


C6 


= 3 


ei 


= 10 


62 


= 20 


63 


= 10 


64 


= 5 


T 


= 10ms 


r 


= 1 




= 0.55 




= 11 




= 0.25 




= 9 


a^cx 


= 1 




= 2 




= 2 




= 2.5 







The strength of connections was scaled by a G [0, 1]. 
Everywhere except in Section IIV CI we used P — 1. In 
Section [IV CI we used either a G [0,1] and /3 = 1, or 
a e [0, 1] and (3 ^ a. 
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